home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / dsyr2k.f < prev    next >
Text File  |  1997-06-25  |  11KB  |  331 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        UPLO, TRANS
  8.       INTEGER            N, K, LDA, LDB, LDC
  9.       DOUBLE PRECISION   ALPHA, BETA
  10. *     .. Array Arguments ..
  11.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  DSYR2K  performs one of the symmetric rank 2k operations
  18. *
  19. *     C := alpha*A*B' + alpha*B*A' + beta*C,
  20. *
  21. *  or
  22. *
  23. *     C := alpha*A'*B + alpha*B'*A + beta*C,
  24. *
  25. *  where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
  26. *  and  A and B  are  n by k  matrices  in the  first  case  and  k by n
  27. *  matrices in the second case.
  28. *
  29. *  Parameters
  30. *  ==========
  31. *
  32. *  UPLO   - CHARACTER*1.
  33. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  34. *           triangular  part  of the  array  C  is to be  referenced  as
  35. *           follows:
  36. *
  37. *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
  38. *                                  is to be referenced.
  39. *
  40. *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
  41. *                                  is to be referenced.
  42. *
  43. *           Unchanged on exit.
  44. *
  45. *  TRANS  - CHARACTER*1.
  46. *           On entry,  TRANS  specifies the operation to be performed as
  47. *           follows:
  48. *
  49. *              TRANS = 'N' or 'n'   C := alpha*A*B' + alpha*B*A' +
  50. *                                        beta*C.
  51. *
  52. *              TRANS = 'T' or 't'   C := alpha*A'*B + alpha*B'*A +
  53. *                                        beta*C.
  54. *
  55. *              TRANS = 'C' or 'c'   C := alpha*A'*B + alpha*B'*A +
  56. *                                        beta*C.
  57. *
  58. *           Unchanged on exit.
  59. *
  60. *  N      - INTEGER.
  61. *           On entry,  N specifies the order of the matrix C.  N must be
  62. *           at least zero.
  63. *           Unchanged on exit.
  64. *
  65. *  K      - INTEGER.
  66. *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
  67. *           of  columns  of the  matrices  A and B,  and on  entry  with
  68. *           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
  69. *           of rows of the matrices  A and B.  K must be at least  zero.
  70. *           Unchanged on exit.
  71. *
  72. *  ALPHA  - DOUBLE PRECISION.
  73. *           On entry, ALPHA specifies the scalar alpha.
  74. *           Unchanged on exit.
  75. *
  76. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
  77. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  78. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  79. *           part of the array  A  must contain the matrix  A,  otherwise
  80. *           the leading  k by n  part of the array  A  must contain  the
  81. *           matrix A.
  82. *           Unchanged on exit.
  83. *
  84. *  LDA    - INTEGER.
  85. *           On entry, LDA specifies the first dimension of A as declared
  86. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  87. *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
  88. *           be at least  max( 1, k ).
  89. *           Unchanged on exit.
  90. *
  91. *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
  92. *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
  93. *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
  94. *           part of the array  B  must contain the matrix  B,  otherwise
  95. *           the leading  k by n  part of the array  B  must contain  the
  96. *           matrix B.
  97. *           Unchanged on exit.
  98. *
  99. *  LDB    - INTEGER.
  100. *           On entry, LDB specifies the first dimension of B as declared
  101. *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
  102. *           then  LDB must be at least  max( 1, n ), otherwise  LDB must
  103. *           be at least  max( 1, k ).
  104. *           Unchanged on exit.
  105. *
  106. *  BETA   - DOUBLE PRECISION.
  107. *           On entry, BETA specifies the scalar beta.
  108. *           Unchanged on exit.
  109. *
  110. *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
  111. *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  112. *           upper triangular part of the array C must contain the upper
  113. *           triangular part  of the  symmetric matrix  and the strictly
  114. *           lower triangular part of C is not referenced.  On exit, the
  115. *           upper triangular part of the array  C is overwritten by the
  116. *           upper triangular part of the updated matrix.
  117. *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  118. *           lower triangular part of the array C must contain the lower
  119. *           triangular part  of the  symmetric matrix  and the strictly
  120. *           upper triangular part of C is not referenced.  On exit, the
  121. *           lower triangular part of the array  C is overwritten by the
  122. *           lower triangular part of the updated matrix.
  123. *
  124. *  LDC    - INTEGER.
  125. *           On entry, LDC specifies the first dimension of C as declared
  126. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  127. *           max( 1, n ).
  128. *           Unchanged on exit.
  129. *
  130. *
  131. *  Level 3 Blas routine.
  132. *
  133. *
  134. *  -- Written on 8-February-1989.
  135. *     Jack Dongarra, Argonne National Laboratory.
  136. *     Iain Duff, AERE Harwell.
  137. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  138. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  139. *
  140. *
  141. *     .. External Functions ..
  142.       LOGICAL            LSAME
  143.       EXTERNAL           LSAME
  144. *     .. External Subroutines ..
  145.       EXTERNAL           XERBLA
  146. *     .. Intrinsic Functions ..
  147.       INTRINSIC          MAX
  148. *     .. Local Scalars ..
  149.       LOGICAL            UPPER
  150.       INTEGER            I, INFO, J, L, NROWA
  151.       DOUBLE PRECISION   TEMP1, TEMP2
  152. *     .. Parameters ..
  153.       DOUBLE PRECISION   ONE         , ZERO
  154.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  155. *     ..
  156. *     .. Executable Statements ..
  157. *
  158. *     Test the input parameters.
  159. *
  160.       IF( LSAME( TRANS, 'N' ) )THEN
  161.          NROWA = N
  162.       ELSE
  163.          NROWA = K
  164.       END IF
  165.       UPPER = LSAME( UPLO, 'U' )
  166. *
  167.       INFO = 0
  168.       IF(      ( .NOT.UPPER               ).AND.
  169.      $         ( .NOT.LSAME( UPLO , 'L' ) )      )THEN
  170.          INFO = 1
  171.       ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
  172.      $         ( .NOT.LSAME( TRANS, 'T' ) ).AND.
  173.      $         ( .NOT.LSAME( TRANS, 'C' ) )      )THEN
  174.          INFO = 2
  175.       ELSE IF( N  .LT.0               )THEN
  176.          INFO = 3
  177.       ELSE IF( K  .LT.0               )THEN
  178.          INFO = 4
  179.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  180.          INFO = 7
  181.       ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN
  182.          INFO = 9
  183.       ELSE IF( LDC.LT.MAX( 1, N     ) )THEN
  184.          INFO = 12
  185.       END IF
  186.       IF( INFO.NE.0 )THEN
  187.          CALL XERBLA( 'DSYR2K', INFO )
  188.          RETURN
  189.       END IF
  190. *
  191. *     Quick return if possible.
  192. *
  193.       IF( ( N.EQ.0 ).OR.
  194.      $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  195.      $   RETURN
  196. *
  197. *     And when  alpha.eq.zero.
  198. *
  199.       IF( ALPHA.EQ.ZERO )THEN
  200.          IF( UPPER )THEN
  201.             IF( BETA.EQ.ZERO )THEN
  202.                DO 20, J = 1, N
  203.                   DO 10, I = 1, J
  204.                      C( I, J ) = ZERO
  205.    10             CONTINUE
  206.    20          CONTINUE
  207.             ELSE
  208.                DO 40, J = 1, N
  209.                   DO 30, I = 1, J
  210.                      C( I, J ) = BETA*C( I, J )
  211.    30             CONTINUE
  212.    40          CONTINUE
  213.             END IF
  214.          ELSE
  215.             IF( BETA.EQ.ZERO )THEN
  216.                DO 60, J = 1, N
  217.                   DO 50, I = J, N
  218.                      C( I, J ) = ZERO
  219.    50             CONTINUE
  220.    60          CONTINUE
  221.             ELSE
  222.                DO 80, J = 1, N
  223.                   DO 70, I = J, N
  224.                      C( I, J ) = BETA*C( I, J )
  225.    70             CONTINUE
  226.    80          CONTINUE
  227.             END IF
  228.          END IF
  229.          RETURN
  230.       END IF
  231. *
  232. *     Start the operations.
  233. *
  234.       IF( LSAME( TRANS, 'N' ) )THEN
  235. *
  236. *        Form  C := alpha*A*B' + alpha*B*A' + C.
  237. *
  238.          IF( UPPER )THEN
  239.             DO 130, J = 1, N
  240.                IF( BETA.EQ.ZERO )THEN
  241.                   DO 90, I = 1, J
  242.                      C( I, J ) = ZERO
  243.    90             CONTINUE
  244.                ELSE IF( BETA.NE.ONE )THEN
  245.                   DO 100, I = 1, J
  246.                      C( I, J ) = BETA*C( I, J )
  247.   100             CONTINUE
  248.                END IF
  249.                DO 120, L = 1, K
  250.                   IF( ( A( J, L ).NE.ZERO ).OR.
  251.      $                ( B( J, L ).NE.ZERO )     )THEN
  252.                      TEMP1 = ALPHA*B( J, L )
  253.                      TEMP2 = ALPHA*A( J, L )
  254.                      DO 110, I = 1, J
  255.                         C( I, J ) = C( I, J ) +
  256.      $                              A( I, L )*TEMP1 + B( I, L )*TEMP2
  257.   110                CONTINUE
  258.                   END IF
  259.   120          CONTINUE
  260.   130       CONTINUE
  261.          ELSE
  262.             DO 180, J = 1, N
  263.                IF( BETA.EQ.ZERO )THEN
  264.                   DO 140, I = J, N
  265.                      C( I, J ) = ZERO
  266.   140             CONTINUE
  267.                ELSE IF( BETA.NE.ONE )THEN
  268.                   DO 150, I = J, N
  269.                      C( I, J ) = BETA*C( I, J )
  270.   150             CONTINUE
  271.                END IF
  272.                DO 170, L = 1, K
  273.                   IF( ( A( J, L ).NE.ZERO ).OR.
  274.      $                ( B( J, L ).NE.ZERO )     )THEN
  275.                      TEMP1 = ALPHA*B( J, L )
  276.                      TEMP2 = ALPHA*A( J, L )
  277.                      DO 160, I = J, N
  278.                         C( I, J ) = C( I, J ) +
  279.      $                              A( I, L )*TEMP1 + B( I, L )*TEMP2
  280.   160                CONTINUE
  281.                   END IF
  282.   170          CONTINUE
  283.   180       CONTINUE
  284.          END IF
  285.       ELSE
  286. *
  287. *        Form  C := alpha*A'*B + alpha*B'*A + C.
  288. *
  289.          IF( UPPER )THEN
  290.             DO 210, J = 1, N
  291.                DO 200, I = 1, J
  292.                   TEMP1 = ZERO
  293.                   TEMP2 = ZERO
  294.                   DO 190, L = 1, K
  295.                      TEMP1 = TEMP1 + A( L, I )*B( L, J )
  296.                      TEMP2 = TEMP2 + B( L, I )*A( L, J )
  297.   190             CONTINUE
  298.                   IF( BETA.EQ.ZERO )THEN
  299.                      C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
  300.                   ELSE
  301.                      C( I, J ) = BETA *C( I, J ) +
  302.      $                           ALPHA*TEMP1 + ALPHA*TEMP2
  303.                   END IF
  304.   200          CONTINUE
  305.   210       CONTINUE
  306.          ELSE
  307.             DO 240, J = 1, N
  308.                DO 230, I = J, N
  309.                   TEMP1 = ZERO
  310.                   TEMP2 = ZERO
  311.                   DO 220, L = 1, K
  312.                      TEMP1 = TEMP1 + A( L, I )*B( L, J )
  313.                      TEMP2 = TEMP2 + B( L, I )*A( L, J )
  314.   220             CONTINUE
  315.                   IF( BETA.EQ.ZERO )THEN
  316.                      C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2
  317.                   ELSE
  318.                      C( I, J ) = BETA *C( I, J ) +
  319.      $                           ALPHA*TEMP1 + ALPHA*TEMP2
  320.                   END IF
  321.   230          CONTINUE
  322.   240       CONTINUE
  323.          END IF
  324.       END IF
  325. *
  326.       RETURN
  327. *
  328. *     End of DSYR2K.
  329. *
  330.       END
  331.